*DECK RINPUT
      SUBROUTINE RINPUT
C==============================================================RINPUT
C      RINPUT READS THE INPUT DATA FILE
C      AND COMPUTES DERIVED QUANTITIES, SCALARS, PLOT DATA,
C      AND THE SPECIES AND REACTION TABULAR DATA.
C----------
C      CALLED BY LAVA
C==============================================================RINPUT
      INCLUDE 'COML.h'
C----------
      INTEGER N,NK,IRE,ICMP,ITMP,IR,ITL,ITR,ITD,ITF,ITB,ITT,NOZ
      INTEGER NPROD,NFCTR
C----------
      DOUBLE PRECISION XMULT,DCHK,HKZERO,TEMPE,
     1     TB,FR,FRL,FRR,FRD,FRF,FRB,FRT
      DOUBLE PRECISION SINTZ,COSTZ,SINTX,COSTX,DOTPRO,ABSPRO
      DOUBLE PRECISION SUMPSZ,RSUMPI
      DOUBLE PRECISION QGAS1,RAD1,RAD2
C
C----------
      EXTERNAL EXITA
      CHARACTER *8 IDSP(NSP)
      CHARACTER *1 IDPH(NSP)
      CHARACTER *8 ID(4)
      CHARACTER *4 MESS
      CHARACTER *50 HDING
      CHARACTER *75 TITLE
C==============================================================RINPUT
C     READ DATA DECK, COMPUTE DERIVED AND SCALAR QUANTITIES
C--------------------------------------------------------------RINPUT
      OPEN(UNIT=NFINP,FILE='lava.inp',STATUS='OLD')
      OPEN(UNIT=NFOUT,FILE='lava.out',STATUS='UNKNOWN')
      REWIND NFINP
      REWIND NFOUT
C--------------------------------------------------------------RINPUT
      READ (NFINP,1020) TITLE
      WRITE(NFOUT,1020) TITLE
      READ (NFINP,1020) TITLE
      WRITE(NFOUT,1020) TITLE
      READ (NFINP,1020) TITLE
      WRITE(NFOUT,1020) TITLE
      READ (NFINP,1020) TITLE
      WRITE(NFOUT,1020) TITLE
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),INCOMP,ID(2),LDT
      WRITE(NFOUT,1100) ID(1),INCOMP,ID(2),LDT
      READ (NFINP,1100) ID(1),ISWIRL,ID(2),ICYL
      WRITE(NFOUT,1100) ID(1),ISWIRL,ID(2),ICYL
C--------------- CONSISTENCY CHECK - GEOMETRIES AND GRIDS -----RINPUT
      IF((NX.EQ.1 .OR. NY.EQ.1 .OR. NZ.GT.1) .AND. ICYL.NE.0)
     1     CALL EXITA(1)
      IF(ICYL.EQ.0 .AND. ISWIRL.NE.0) CALL EXITA(1)
C--------------------------------------------------------------RINPUT
      READ (NFINP,1100) ID(1),NODIFF,ID(2),LWALL
      WRITE(NFOUT,1100) ID(1),NODIFF,ID(2),LWALL
      READ (NFINP,1100) ID(1),NLTE
      WRITE(NFOUT,1100) ID(1),NLTE
      READ (NFINP,1100) ID(1),ICONDP
      WRITE(NFOUT,1100) ID(1),ICONDP
      READ (NFINP,1100) ID(1),IEVAP
      WRITE(NFOUT,1100) ID(1),IEVAP
      READ (NFINP,1100) ID(1),IPOST
      WRITE(NFOUT,1100) ID(1),IPOST
      READ (NFINP,1100) ID(1),IFLUCT
      WRITE(NFOUT,1100) ID(1),IFLUCT
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1110) ID(1),NMAX,ID(2),NCYC,ID(3),JUMP
      WRITE(NFOUT,1110) ID(1),NMAX,ID(2),NCYC,ID(3),JUMP
      READ (NFINP,1120) ID(1),DTJUMP
      WRITE(NFOUT,1130) ID(1),DTJUMP
      READ (NFINP,1120) ID(1),DTI,ID(2),DTMAX,ID(3),TMAX
      WRITE(NFOUT,1130) ID(1),DTI,ID(2),DTMAX,ID(3),TMAX
      READ (NFINP,1120) ID(1),TIMPTL
      WRITE(NFOUT,1130) ID(1),TIMPTL
C ----- input of the safty factor for timestep by Y.P. WAN
      READ (NFINP,1120) ID(1),SAFEDT,ID(2),SAFDTD,ID(3),DTGROW
      WRITE(NFOUT,1130) ID(1),SAFEDT,ID(2),SAFDTD,ID(3),DTGROW
C ------------------------------- GRID CONTROL by Y.P. WAN
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1110) ID(1),NX,ID(2),NY,ID(3),NZ
      WRITE(NFOUT,1110) ID(1),NX,ID(2),NY,ID(3),NZ
      READ (NFINP,1120) ID(1),X(1),ID(2),Y(1),ID(3),Z(1)
      WRITE(NFOUT,1130) ID(1),X(1),ID(2),Y(1),ID(3),Z(1)
      READ (NFINP,1120) ID(1),XL,ID(2),YL,ID(3),ZL
      WRITE(NFOUT,1130) ID(1),XL,ID(2),YL,ID(3),ZL
      READ (NFINP,1110) ID(1),IGRIDX,ID(2),IGRIDY,ID(3),IGRIDZ
      WRITE(NFOUT,1110) ID(1),IGRIDX,ID(2),IGRIDY,ID(3),IGRIDZ
      READ (NFINP,1110) ID(1),ICIN,ID(2),ICOUT
      WRITE(NFOUT,1110) ID(1),ICIN,ID(2),ICOUT
C------------------------- CONSISTENCY CHECK - DIMENSIONS -----RINPUT
      IF(NXT.EQ.1 .OR. NX.EQ.1) THEN
        IF(NXT.NE.1 .OR. NX.NE.1) CALL EXITA(1)
        IF(XL.LE.ZERO) CALL EXITA(1)
        IF(IGRIDX.NE.0) CALL EXITA(1)
      ELSE
        IF(NXT.NE.NX+2) CALL EXITA(1)
      END IF
      IF(NYT.EQ.1 .OR. NY.EQ.1) THEN
        IF(NYT.NE.1 .OR. NY.NE.1) CALL EXITA(1)
        IF(YL.LE.ZERO) CALL EXITA(1)
        IF(IGRIDY.NE.0) CALL EXITA(1)
      ELSE
        IF(NYT.NE.NY+2) CALL EXITA(1)
      END IF
      
      IF(NZT.EQ.1 .OR. NZ.EQ.1) THEN
        IF(NZT.NE.1 .OR. NZ.NE.1) CALL EXITA(1)
        IF(ZL.LE.ZERO) CALL EXITA(1)
        IF(IGRIDZ.NE.0) CALL EXITA(1)
      ELSE
        IF(NZT.NE.NZ+2) CALL EXITA(1)
      END IF
C------------------------------------------- ACCELERATION -----RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),MACCEL
      WRITE(NFOUT,1100) ID(1),MACCEL
      READ (NFINP,1120) ID(1),PGS
      WRITE(NFOUT,1130) ID(1),PGS
           RPGS2=ONE/(PGS*PGS)
      READ (NFINP,1120) ID(1),DAMP,ID(2),THETA,ID(3),RELAX,ID(4),PQDIF
      WRITE(NFOUT,1130) ID(1),DAMP,ID(2),THETA,ID(3),RELAX,ID(4),PQDIF
      PREF=ZERO
      IF(MACCEL.EQ.1) PREF=PAMB
C--------------- CONSISTENCY CHECK - ACCELERATION METHODS -----RINPUT
      IF(MACCEL.LT.0 .OR. MACCEL.GT.2) CALL EXITA(1)
      IF(MACCEL.NE.0 .AND. PGS.NE.ONE) CALL EXITA(1)
      IF(MACCEL.EQ.0 .AND. LDT.EQ.2) CALL EXITA(1)
c--- till we get relation for dive()
      if(maccel.eq.0 .and. damp.ne.zero) CALL EXITA(1)
      IF(THETA.GE.ONE) CALL EXITA(1)
      IF(ICYL.EQ.1  .AND. GZ.NE.ZERO) CALL EXITA(1)
      DCHK=HALF*(ONE-THETA)
      IF(MACCEL.NE.0 .AND. DAMP.GE.DCHK) CALL EXITA(1)
      IF(INCOMP.EQ.1) THEN
        IF(MACCEL.NE.2 .OR. NSP.NE.1) CALL EXITA(1)
      END IF
      IF(MACCEL.EQ.2 .AND. NLTE.NE.0) CALL EXITA(1)
C--------------------------------------------- DONOR-CELL -----RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1120) ID(1),ADC,ID(2),BDC
      WRITE(NFOUT,1130) ID(1),ADC,ID(2),BDC
C-------------------------------------- TURBULENCE MODELS -----RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),ITURB
      WRITE(NFOUT,1100) ID(1),ITURB
      READ (NFINP,1100) ID(1),KECOR
      WRITE(NFOUT,1100) ID(1),KECOR
C--------------- CONSISTENCY CHECK - K-EPSILON CORRECTION -----RINPUT
      IF(ICYL.NE.1 .AND. KECOR.EQ.1) CALL EXITA(1)
C==============================================================RINPUT
C
C     TORCH OPERATING CONDITIONS, ADDED BY Y.P. WAN 10-19-97
C
C==============================================================RINPUT
C     CONDITIONS ARE GIVEN FROM THE INPUT FILE INSTEAD OF IN START
C     ALL THE VARIABLES INTRODUCED HERE SHOULD BE INCLUDED IN COML.H
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1120) ID(1),CURCY
      WRITE(NFOUT,1130) ID(1),CURCY
      READ (NFINP,1120) ID(1),VOLTG
      WRITE(NFOUT,1130) ID(1),VOLTG
      READ (NFINP,1120) ID(1),EFFCY
      WRITE(NFOUT,1130) ID(1),EFFCY
      READ (NFINP,1120) ID(1),FLRT1
      WRITE(NFOUT,1130) ID(1),FLRT1
      READ (NFINP,1120) ID(1),FLRT2
      WRITE(NFOUT,1130) ID(1),FLRT2
      READ (NFINP,1120) ID(1),PWVEL
      WRITE(NFOUT,1130) ID(1),PWVEL
      READ (NFINP,1120) ID(1),PWTEP
      WRITE(NFOUT,1130) ID(1),PWTEP
      READ (NFINP,1120) ID(1),SWIRN
      WRITE(NFOUT,1130) ID(1),SWIRN
C================== AMBIENT CONDITIONS BY Y.P. WAN 10-21-97
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1120) ID(1),PAMB
      WRITE(NFOUT,1130) ID(1),PAMB
      READ (NFINP,1120) ID(1),TEMAMB
      WRITE(NFOUT,1130) ID(1),TEMAMB
      READ (NFINP,1120) ID(1),GX
      WRITE(NFOUT,1130) ID(1),GX
      READ (NFINP,1120) ID(1),GY
      WRITE(NFOUT,1130) ID(1),GY
      READ (NFINP,1120) ID(1),GZ
      WRITE(NFOUT,1130) ID(1),GZ
C==============================================================RINPUT
C
C     SPECIES
C
C==============================================================RINPUT
C     THE FORMULA, INITIAL DENSITY, MOLECULAR WEIGHT, AND HEAT
C     OF FORMATION ARE READ IN BELOW.  HTFORM IS CONVERTED
C     FROM KCAL/MOLE TO CGS UNITS FOR ALL THE SPECIES.
C     NSP SHOULD BE GIVEN AS A PARAMETER
C--------------------- NSP SHOULD BE GIVEN AS A PARAMETER -----RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
C---------- INDEX FOR ELECTRON (IELC = 0 FOR NO ELECTRON)
      READ (NFINP,1100) ID(1),IELC
      WRITE(NFOUT,1100) ID(1),IELC
C---------------------- CONSISTENCY CHECK - NON-LTE MODEL -----RINPUT
      IF(NLTE.NE.0 .AND. IELC.EQ.0) CALL EXITA(1)
C--------------------------------------------------------------RINPUT
      DO 200 ISP=1,NSP
      READ (NFINP,1160) IDSP(ISP),IDPH(ISP),ID(1),MW(ISP),
     1                  ID(2),HTFORM(ISP),ID(3),QOM(ISP)
      IF(IDPH(ISP).EQ.'G') THEN
        IGAS(ISP)=NGAS
      ELSEIF(IDPH(ISP).EQ.'C') THEN
        IGAS(ISP)=NCNDNS
      ELSE
        CALL EXITA(1)
      END IF
      RMW(ISP)=ONE/MW(ISP)
      HTFORM(ISP)=HTFORM(ISP)*ERGCAL
      QOM(ISP)=QOM(ISP)*ELECHG*AVOGAD*RMW(ISP)
      WRITE(NFOUT,1170) IDSP(ISP),IDPH(ISP),ID(1),MW(ISP),
     1                  ID(2),HTFORM(ISP),ID(3),QOM(ISP)
  200 CONTINUE
C==============================================================RINPUT
C
C     CHEMISTRY LOGIC
C
C========================== READ IN KINETIC REACTION DATA =====RINPUT
C     REACTIONS INVOLVING MULTIPLE THIRD BODY CAN BE TREATED AS A
C     SINGLE REACTION USING THIRD BODY CHAPERONE EFFICIENCIES.
C     WHEN CHAPERONE EFFICIENCY IS USED FOR THIRD BODY, DO NOT PUT
C     THIRD BODY IN STOICHIOMETIRC COEFFICIENTS.
C--------------------------------------------------------------RINPUT
C     WHEN THERE IS FORWARD REACTION ONLY, PUT LARGE NUMBER IN CKS.
C     SWITCH FORWARD AND BACKWARD FOR BACKWARD REACTION ONLY CASE.
C--------------------------------------------------------------RINPUT
C     TFLAGK=1.0 FOR REACTIONS DEPEND ON ELECTRON TEMPERTURE
C            0.0 FOR REACTIONS DEPEND ON HEAVY PARTICLE TEMPERTURE
C     FRK=FRACTION OF TOTAL CHEMICAL ENERGY RELEASE LOST AS RADIATION
C     FEK=FRACTION OF TOTAL CHEMICAL ENERGY RELEASE THAT GOES TO
C         ELECTRON
C     FEK + FRK <= 1.0
C     FEK + FRK = 1.0 FOR TFLAGK = 1.0
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),NRK
      WRITE(NFOUT,1100) ID(1),NRK
      IF(NRK.EQ.0) GO TO 300
      DO 220 IR=1,NRK
      READ (NFINP,1010) HDING
      READ (NFINP,1210) FEK(IR),FRK(IR),TFLAGK(IR)
      IF((FEK(IR)+FRK(IR)).GT.ONE) CALL EXITA(1)
      IF(TFLAGK(IR).EQ.ONE .AND. (FEK(IR)+FRK(IR)).LT.(ONE-SMALL))
     1      CALL EXITA(1)
      READ (NFINP,1210) CF(IR),EF(IR),ZETAF(IR)
      IF(CF(IR).EQ.ZERO) CALL EXITA(1)
      READ (NFINP,1220) AKS(IR),BKS(IR),CKS(IR),DKS(IR),EKS(IR)
      READ (NFINP,1230) ID(1),TCUTL(IR),ID(2),TCUTH(IR)
      READ (NFINP,1240) (AM(ISP,IR),ISP=1,NSP)
      READ (NFINP,1240) (BM(ISP,IR),ISP=1,NSP)
      READ (NFINP,1250) (CHAEFF(ISP,IR),ISP=1,NSP)
      READ (NFINP,1250) (AE(ISP,IR),ISP=1,NSP)
      READ (NFINP,1250) (BE(ISP,IR),ISP=1,NSP)
C-------------------------- REACTION-SPECIES INDEX ARRAYS -----RINPUT
      NK=0
      QR(IR)=ZERO
      DO 210 ISP=1,NSP
      IF(AM(ISP,IR).EQ.0 .AND. BM(ISP,IR).EQ.0) GO TO 210
      NK=NK+1
      CM(NK,IR)=ISP
      FAM(ISP,IR)=DBLE(AM(ISP,IR))
      FBM(ISP,IR)=DBLE(BM(ISP,IR))
      FBMAM(ISP,IR)=FBM(ISP,IR)-FAM(ISP,IR)
      QR(IR)=QR(IR)-FBMAM(ISP,IR)*HTFORM(ISP)
  210 CONTINUE
      NELEM(IR)=NK
      WRITE(NFOUT,1010) HDING
      WRITE(NFOUT,1300) IR,QR(IR)
      WRITE(NFOUT,1310) FEK(IR),FRK(IR),TFLAGK(IR)
      WRITE(NFOUT,1320) CF(IR),EF(IR),ZETAF(IR)
      WRITE(NFOUT,1330) AKS(IR),BKS(IR),CKS(IR),DKS(IR),EKS(IR)
      WRITE(NFOUT,1340) ID(1),TCUTL(IR),ID(2),TCUTH(IR)
      MESS=' LHS'
      WRITE(NFOUT,1350) MESS,(AM(ISP,IR),ISP=1,NSP)
      MESS=' RHS'
      WRITE(NFOUT,1350) MESS,(BM(ISP,IR),ISP=1,NSP)
      WRITE(NFOUT,1360) (CHAEFF(ISP,IR),ISP=1,NSP)
      MESS=' FOR'
      WRITE(NFOUT,1370) MESS,(AE(ISP,IR),ISP=1,NSP)
      MESS='BACK'
      WRITE(NFOUT,1370) MESS,(BE(ISP,IR),ISP=1,NSP)
C---------------------- CONSISTENCY CHECK - NON-LTE MODEL -----RINPUT
      IF(NLTE.EQ.0 .AND. FEK(IR).NE.ZERO) CALL EXITA(1)
      IF(NLTE.EQ.0 .AND. TFLAGK(IR).NE.ZERO) CALL EXITA(1)
  220 CONTINUE
C===== READ IN FAST KINETIC AND EQUILIBRIUM REACTION DATA =====RINPUT
C     WHEN REACTIONS NEED TO BE TREATED AS EQUILIBRIUM, PUT LARGE
C     NUMBER IN CFE.  FAST KINETIC REACTIONS ARE ASSUMED TO BE
C     ELEMENTARY REACTIONS (STOICHIOMETIRC COEFFIECIENTS USED FOR
C     EXPONENTS) AND ALWAYS HAVE FORWARD AND BACKWARD REACTIONS.
C--------------------------------------------------------------RINPUT
C     TFLAGE=1.0 FOR REACTIONS DEPEND ON ELECTRON TEMPERTURE
C            0.0 FOR REACTIONS DEPEND ON HEAVY PARTICLE TEMPERTURE
C     FRE=FRACTION OF TOTAL CHEMICAL ENERGY RELEASE LOST AS RADIATION
C     FEE=FRACTION OF TOTAL CHEMICAL ENERGY RELEASE THAT GOES TO
C         ELECTRON
C     FEE + FRE <= 1.0
C     FEE + FRE = 1.0 FOR TFLAGE = 1.0
C--------------------------------------------------------------RINPUT
  300 READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),NRE
      WRITE(NFOUT,1100) ID(1),NRE
      IF(INCOMP.EQ.1) THEN
        IF(NRK.NE.0 .OR. NRE.NE.0) CALL EXITA(1)
      END IF
      IF(NRE.EQ.0) GO TO 400
      DO 330 IRE=1,NRE
      READ (NFINP,1010) HDING
      READ (NFINP,1210) FEE(IRE),FRE(IRE),TFLAGE(IRE)
      IF((FEE(IRE)+FRE(IRE)).GT.ONE) CALL EXITA(1)
      IF(TFLAGE(IRE).EQ.ONE .AND. (FEE(IRE)+FRE(IRE)).LT.(ONE-SMALL))
     1      CALL EXITA(1)
      READ (NFINP,1210) CFE(IRE),EFE(IRE),ZETAFE(IRE)
      IF(CFE(IRE).EQ.ZERO) CALL EXITA(1)
      READ (NFINP,1220) AS(IRE),BS(IRE),CS(IRE),DS(IRE),ES(IRE)
      READ (NFINP,1230) ID(1),TCUTEL(IRE),ID(2),TCUTEH(IRE)
      READ (NFINP,1240) (AN(ISP,IRE),ISP=1,NSP)
      READ (NFINP,1240) (BN(ISP,IRE),ISP=1,NSP)
C--------------------------------------------------------------RINPUT
      NK=0
      QEQ(IRE)=ZERO
      DO 310 ISP=1,NSP
        FBNAN(ISP,IRE)=DBLE(BN(ISP,IRE)-AN(ISP,IRE))
        IF(AN(ISP,IRE).EQ.0 .AND. BN(ISP,IRE).EQ.0) GO TO 310
        NK=NK+1
        CN(NK,IRE)=ISP
        QEQ(IRE)=QEQ(IRE)-FBNAN(ISP,IRE)*HTFORM(ISP)
  310 CONTINUE
      NLM(IRE)=NK
      WRITE(NFOUT,1405) FEE(IRE),FRE(IRE),TFLAGE(IRE)
      WRITE(NFOUT,1400) CFE(IRE),EFE(IRE),ZETAFE(IRE)
      WRITE(NFOUT,1410)
     1     IRE,HDING,AS(IRE),BS(IRE),CS(IRE),DS(IRE),ES(IRE)
      WRITE(NFOUT,1340) ID(1),TCUTEL(IRE),ID(2),TCUTEH(IRE)
      WRITE(NFOUT,1420) QEQ(IRE),NLM(IRE)
      MESS='AN'
      WRITE(NFOUT,1430) MESS,(AN(ISP,IRE),ISP=1,NSP)
      MESS='BN'
      WRITE(NFOUT,1430) MESS,(BN(ISP,IRE),ISP=1,NSP)
      MESS='CN'
      WRITE(NFOUT,1430) MESS,(CN(ISP,IRE),ISP=1,NSP)
  330 CONTINUE
C==============================================================RINPUT
C     CONVERT HK ARRAYS FROM ENTHALPY TO EK ARRAYS OF SIE (ERGS/GRAM)
C--------------------------------------------------------------RINPUT
  400 DO 410 ISP=1,NSP
      HKZERO=HK(1,ISP)
      DO 410 N=1,NIT
      TEMPE=HUNDRD*DBLE(N-1)
      EK(N,ISP)=RMW(ISP)*(ERGCAL*(HK(N,ISP)-HKZERO)-RGAS*TEMPE)
  410 CONTINUE
C==============================================================RINPUT
C
C     PARTICLE LOGIC
C
C==============================================================RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),IPTCL,ID(2),I3DP2D,ID(3),ICARRY
      WRITE(NFOUT,1100) ID(1),IPTCL,ID(2),I3DP2D,ID(3),ICARRY
C--------------------------------------------------------------RINPUT
      IF(IPTCL.EQ.0) GO TO 1000
C--------------------------------------------------------------RINPUT
C     DO NOT USE PARTICLE ON 1-D
C--------------------------------------------------------------RINPUT
      IF(NX.EQ.1 .AND. IPTCL.EQ.1) CALL EXITA(1)
      IF(NY.EQ.1 .AND. IPTCL.EQ.1) CALL EXITA(1)
C--------------------------------------------------------------RINPUT
C     I3DP2D = 1 (3-D PARTICLE IN 2-D AXISYMMETRIC FLOW) SHOULD NOT
C     BE USED FOR 3-D, SINCE I3DP2D IS USED IN COORDINATE
C     TRANSFORMATION.
C--------------------------------------------------------------RINPUT
C     POSSIBLE CASES
C   NZ=1, ICYL=0, ISWIRL=0, I3DP2D=0: PLANE 2-D
C   NZ=1, ICYL=1, ISWIRL=0, I3DP2D=0: CYL, NO SWIRL, 2-D PARTICLE
C   NZ=1, ICYL=1, ISWIRL=1, I3DP2D=0: CYL, SWIRL, 2-D PARTICLE
C          -- DO NOT EXIST!!
C   NZ=1, ICYL=1, ISWIRL=0, I3DP2D=1: CYL, NO SWIRL, 3-D PART, ZINJ=0
C   NZ=1, ICYL=1, ISWIRL=1, I3DP2D=1: CYL, SWIRL, 3-D PARTICLE
C   NZ>1, ICYL=0, ISWIRL=0, I3DP2D=0: PLANE 3-D
C-------------------------------------- CONSISTENCY CHECK -----RINPUT
      IF(LDT.EQ.2 .AND. IPTCL.EQ.1) CALL EXITA(1)
      IF(I3DP2D.EQ.1 .AND. ICYL.NE.1) CALL EXITA(1)
      IF(ISWIRL.EQ.1 .AND. I3DP2D.EQ.0) CALL EXITA(1)
C==============================================================RINPUT
C     READ IN PARTICLE MATERIAL PROPERTIES
C--------------------------------------------------------------RINPUT
C     ROPI = PARTICLE DENSITY (GRAM/CM**3)
C     TMMI = MINIMUM MELTING TEMPERATURE
C     TMLI = MAXIMUM MELTING TEMPERATURE, TML = TMM FOR PURE SUBSTANCE
C     HLTNT = LATENT HEAT OF FUSION (CAL/G -> ERG/G)
C           = EP(TMLI) - EP(TMMI)
C             (FOR MIXTURES, T CHANGE DURING MELTING)
C     PSPHS = SPECIFIC HEAT (SOLID PHASE: CAL/(GRAM*K) -> ERG/(GRAM*K))
C     PSPHL = SPECIFIC HEAT (LIQUID PHASE: CAL/(GRAM*K) -> ERG/(GRAM*K))
C     EMSSPI = TOTAL EMISSIVITY (ZERO TO NEGLECT RADIATION FROM PARTICLE)
C--------------------------------------------------------------RINPUT
C     NOTE THAT RPSPHS AND RPSPHL ARE USED FOR LINEAR EQUATION OF
C     STATE OF PARTICLES. IF OTHER MODELS ARE USED, THIS PART MAY
C     HAS TO BE MODIFIED (SUBROUTINE PSTATE SHOULD BE ALSO ALTERED
C     ACCORDINGLY).
C--------------------------------------------------------------RINPUT
      DO 540 IPSP=1,NPSP
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1220) ROPI(IPSP),TMMI(IPSP),TMLI(IPSP),HLTNT(IPSP)
      WRITE(NFOUT,1600) ROPI(IPSP),TMMI(IPSP),TMLI(IPSP),HLTNT(IPSP)
      READ (NFINP,1220) PSPHS(IPSP),PSPHL(IPSP),EMSSPI(IPSP)
      WRITE(NFOUT,1610) PSPHS(IPSP),PSPHL(IPSP),EMSSPI(IPSP)
C--------------------------------------------------------------RINPUT
C     CONVERT HLTNT, PSHPHS, PSPHL TO ERG/GRAM, ERG/GRAM/K, ERG/GRAM/K
C--------------------------------------------------------------RINPUT
      HLTNT(IPSP)=HLTNT(IPSP)*ERGCAL*THOUTH
      PSPHS(IPSP)=PSPHS(IPSP)*ERGCAL*THOUTH
      PSPHL(IPSP)=PSPHL(IPSP)*ERGCAL*THOUTH
  540 CONTINUE
C-------------------------------------- NOZZLE INJECTIONS -----RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1100) ID(1),NUMNOZ
      WRITE(NFOUT,1100) ID(1),NUMNOZ
      IF(NUMNOZ.EQ.0) GO TO 1000
C--------------------------------------------------------------RINPUT
      DO 530 NOZ=1,NUMNOZ
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
C--------------------------------------------------------------RINPUT
C     ZINJ IN 2-D IS TRUE Z POSITION, NOT R*THETA
C     ZINJ IN 2-D DOES NOT HAVE TO BE ZERO.  OFF AXIS PLANE INJECTION
C     IS ALLOWD IN 3-D PARTICLE IN 2-D FLOW WITH SWIRL SITUATIONS.
C--------------------------------------------------------------RINPUT
C     TSPMAS = MASS INJECTION RATE (GRAM/SEC)
C     ANGLE --- DEGREE
C--------------------------------------------------------------RINPUT
C     FOR 2-D, ANOZ IS (UNIT DEPTH)*(1-D NOZZLE SIZE), AND LOADING
C     RATES FOR PARTICLES AND GASES ARE PER UNIT DEPTH (CM).
C     FOR 2-D, AXISYMMETIRC CASE, TOTAL LOADING RATE, AND TOTAL
C     NOZZLE AREA (CIRCUMFERENTIAL).
C     FOR 2-D, NO SWIRL, TILTZ = 90 DEGREE (PLANE AND AXISYMMETRIC)
C--------------------------------------------------------------RINPUT
      READ (NFINP,1220) XINJ(NOZ),YINJ(NOZ),ZINJ(NOZ)
      WRITE(NFOUT,1500) XINJ(NOZ),YINJ(NOZ),ZINJ(NOZ)
      READ (NFINP,1220) CONE(NOZ),TILTZ(NOZ),TILTX(NOZ)
      WRITE(NFOUT,1530) CONE(NOZ),TILTZ(NOZ),TILTX(NOZ)
      READ (NFINP,1510) ANOZ(NOZ),T1INJ(NOZ)
      WRITE(NFOUT,1520) ANOZ(NOZ),T1INJ(NOZ)
      READ (NFINP,1510) TSPMAS(NOZ),PMINJ(NOZ)
      WRITE(NFOUT,1540) TSPMAS(NOZ),PMINJ(NOZ)
C--------------------------------------------------------------RINPUT
C     INJECTION VELOCITY IS AN INPUT VARIABLE: PARTICLES ARE NOT
C     FULLY RELAXED WITH CARRIER GAS IN THE SMALL INJECTION TUBE,
C     J. R. FINCKE, PRIVATE COMMUNICATION (1996).
C--------------------------------------------------------------RINPUT
      READ (NFINP,1220) VELINJ(NOZ),VDISP(NOZ),TPINJ(NOZ)
      WRITE(NFOUT,1550) VELINJ(NOZ),VDISP(NOZ),TPINJ(NOZ)
C--------------------------------------------------------------RINPUT
C     PYINJ = PARTICLE MASS FRACTION OF SPECIES IPSP
C     PXINJ = PARTICLE VOLUME FRACTION OF SPECIES IPSP
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1250) (PYINJ(NOZ,IPSP),IPSP=1,NPSP)
      WRITE(NFOUT,1605) (PYINJ(NOZ,IPSP),IPSP=1,NPSP)
C---------------------------------------- CALCULATE PXINJ -----RINPUT
      RSUMPI=ZERO
      DO 550 IPSP=1,NPSP
        RSUMPI=RSUMPI+PYINJ(NOZ,IPSP)/ROPI(IPSP)
  550 CONTINUE
      RSUMPI=ONE/RSUMPI
      DO 560 IPSP=1,NPSP
        PXINJ(NOZ,IPSP)=RSUMPI*PYINJ(NOZ,IPSP)/ROPI(IPSP)
  560 CONTINUE
C--------------------------------------------------------------RINPUT
C     RPMIN = MINIMUM PARTICLE RADIUS (CM)
C     RPMAX = MAXIMUM PARTICLE RADIUS (CM)
C     NPSDIS = NUMBER OF INTERVALS ON PARTICLE SIZE DISTRIBUTION
C     PSZDIS = PARTICLE SIZE DISTRIBUTION - NEED NOT BE NORMALIZED
C              NEED TO HAVE UNIFORM INTERVAL
C     PSZDEL = PARTICLE SIZE DISTRIBUTION INTERVAL
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      READ (NFINP,1220) RPMIN(NOZ),RPMAX(NOZ)
      WRITE(NFOUT,1615) RPMIN(NOZ),RPMAX(NOZ)
      READ (NFINP,1100) ID(1),NPSDIS(NOZ)
      WRITE(NFOUT,1100) ID(1),NPSDIS(NOZ)
      IF(NPSDIS(NOZ).GT.NSIZE) CALL EXITA(1)
      IF(NPSDIS(NOZ).LE.5) CALL EXITA(1)
      PSZDEL(NOZ)=(RPMAX(NOZ)-RPMIN(NOZ))/DBLE(NPSDIS(NOZ))
C--------------------------------------------------------------RINPUT
C     PSZINT = INTEGRAL OF PSZDIS*RADP**3 - NORMALIZED
C     IF EFFECTIVE RADIUS IS USED, PARTICLES NEED NOT BE SPHERICAL
C          (SEE INJECT)
C--------------------------------------------------------------RINPUT
      PSZINT(1,NOZ)=ZERO
      RAD1=RPMIN(NOZ)
        DO 510 I=1,NPSDIS(NOZ)
        READ (NFINP,1560) PSZDIS(I,NOZ)
        RAD2=RAD1+PSZDEL(NOZ)
        PSZINT(I+1,NOZ)=PSZINT(I,NOZ)+PSZDIS(I,NOZ)*(RAD2**4-RAD1**4)
        RAD1=RAD2
  510   CONTINUE
      SUMPSZ=PSZINT(NPSDIS(NOZ)+1,NOZ)
        DO 520 I=1,NPSDIS(NOZ)+1
        PSZINT(I,NOZ)=PSZINT(I,NOZ)/SUMPSZ
        WRITE(NFOUT,1570) RPMIN(NOZ),RPMIN(NOZ)+DBLE(I-1)*PSZDEL(NOZ),
     1                    PSZINT(I,NOZ)
  520   CONTINUE
C--------------------------- UNIT VECTORS FOR EACH NOZZLE -----RINPUT
C     CONVERT CONE, TILTZ, AND TILTX FROM DEGREE TO RADIAN
C--------------------------------------------------------------RINPUT
      CONE(NOZ)=CONE(NOZ)*PIO180
      TILTZ(NOZ)=TILTZ(NOZ)*PIO180
      TILTX(NOZ)=TILTX(NOZ)*PIO180
      COSTZ=COS(TILTZ(NOZ))
      SINTZ=SIN(TILTZ(NOZ))
      COSTX=COS(TILTX(NOZ))
      SINTX=SIN(TILTX(NOZ))
      EAVEC(NOZ,1)=SINTZ*COSTX
      EAVEC(NOZ,2)=SINTZ*SINTX
      EAVEC(NOZ,3)=COSTZ
      ENVEC(NOZ,1)=SIN(TILTZ(NOZ)+PI*HALF)*COSTX
      ENVEC(NOZ,2)=SIN(TILTZ(NOZ)+PI*HALF)*SINTX
      ENVEC(NOZ,3)=COS(TILTZ(NOZ)+PI*HALF)
      EOVEC(NOZ,1)=EAVEC(NOZ,2)*ENVEC(NOZ,3)-EAVEC(NOZ,3)*ENVEC(NOZ,2)
      EOVEC(NOZ,2)=EAVEC(NOZ,3)*ENVEC(NOZ,1)-EAVEC(NOZ,1)*ENVEC(NOZ,3)
      EOVEC(NOZ,3)=EAVEC(NOZ,1)*ENVEC(NOZ,2)-EAVEC(NOZ,2)*ENVEC(NOZ,1)
C------------------ CONSISTENCY CHECK - NOZZLE INJECTIONS -----RINPUT
C     ONLY INJECTION TOWARD THE AXIS IS ALLOWD IN 3-D PARTICLE IN 2-D
C     FLOW WITHOUT SWIRL
C--------------------------------------------------------------RINPUT
      IF(ANOZ(NOZ).LE.ZERO) CALL EXITA(1)
      IF(ISWIRL.EQ.0 .AND. NZ.EQ.1) THEN
        DOTPRO=ABS(EAVEC(NOZ,1)*XINJ(NOZ)+EAVEC(NOZ,3)*ZINJ(NOZ))
        ABSPRO=SQRT(XINJ(NOZ)**2+ZINJ(NOZ)**2)
     1       *SQRT(EAVEC(NOZ,1)**2+EAVEC(NOZ,3)**2)
        IF(ONE-(DOTPRO/ABSPRO).GT.SMALL) CALL EXITA(1)
      END IF
C--------------------------------------------------------------RINPUT
  530 CONTINUE
C==============================================================RINPUT
C
C     CARRIER GAS LOGIC
C
C==============================================================RINPUT
      IF(ICARRY.EQ.0) GO TO 1000
C--------------------------------------------------------------RINPUT
C     CARRYG - GRAM/SEC
C     WE NEED TO MODIFY BELOW TO HAVE NON-LTE OPTION FOR CARRIER GAS.
C--------------------------------------------------------------RINPUT
      READ (NFINP,1010) HDING
      WRITE(NFOUT,1010) HDING
      DO 700 ISP=1,NSP
      READ (NFINP,1220) (CARRYG(NOZ,ISP),NOZ=1,NUMNOZ)
      WRITE(NFOUT,1620) ISP,(CARRYG(NOZ,ISP),NOZ=1,NUMNOZ)
C------------ CONSISTENCY CHECK - CARRIER GAS AND NON-LTE -----RINPUT
      IF(NLTE.NE.0 .AND. CARRYG(NOZ,IELC).NE.ZERO) CALL EXITA(1)
  700 CONTINUE
C--------------------------------------------------------------RINPUT
C     MOMENTUM SOURCES (CRMOMX, CRMOMY, AND CRMOMZ) ARE STORED
C     CRMOMX,Y,Z = (MOMENTUM SOURCE) * PAMB
C     THESE ARE CORRECTED AT SUBROUTINE CARRY USING LOCAL PRESSURE
C     ENERGY SOURCE IS STORED AT CRENGY
C--------------------------------------------------------------RINPUT
C     TEMPERATURE ARE ASSUMED TO BE THE SAME
C     FOR CARRIER GASES AND PARTICLES FOR EACH NOZZLE
C     TEMPERATURE AND PRESSURE ARE THE SAME AS AMBIENT CONDITION
C--------------------------------------------------------------RINPUT
      DO 720 NOZ=1,NUMNOZ
      QGAS1=ZERO
      CRGAS(NOZ)=ZERO
      CRENGY(NOZ)=ZERO
      TB=HUNDTH*TPINJ(NOZ)
      IT=INT(TB)
      FR=TB-DBLE(IT)
        DO 710 ISP=1,NSP
        QGAS1=QGAS1+CARRYG(NOZ,ISP)*RMW(ISP)
        CRGAS(NOZ)=CRGAS(NOZ)+CARRYG(NOZ,ISP)
        CRENGY(NOZ)=CRENGY(NOZ)+CARRYG(NOZ,ISP)
     1              *((ONE-FR)*EK(IT+1,ISP)+FR*EK(IT+2,ISP))
  710   CONTINUE
      CRMOMX(NOZ)=HALF*CRGAS(NOZ)*VELINJ(NOZ)*EAVEC(NOZ,1)*PAMB
      CRMOMY(NOZ)=HALF*CRGAS(NOZ)*VELINJ(NOZ)*EAVEC(NOZ,2)*PAMB
      CRMOMZ(NOZ)=HALF*CRGAS(NOZ)*VELINJ(NOZ)*EAVEC(NOZ,3)*PAMB
C--------------------------------------------------------------RINPUT
  720 CONTINUE
C--------------------------------------------------------------RINPUT
 1000 CONTINUE
C==============================================================RINPUT
      RETURN
 1010 FORMAT(1X,A50)
 1020 FORMAT(A75)
 1100 FORMAT(A8,I5)
 1110 FORMAT(A8,I7)
 1120 FORMAT(A8,F10.5)
 1130 FORMAT(A8,2X,1PE12.5)
 1160 FORMAT(A8,2X,A1,4(A6,F10.4))
 1170 FORMAT(A8,2X,A1,4(A6,1PE14.5))
 1210 FORMAT(3(6X,F10.5))
 1220 FORMAT(5(6X,F10.5))
 1230 FORMAT(2(A9,F10.5))
 1240 FORMAT(5X,15I5)
 1250 FORMAT(5X,8F8.3)
 1300 FORMAT(//32X,'KINETIC REACTION',I5/' QR =' ,1PE20.6)
 1310 FORMAT(' FEK, FRK, TFLAGK =',1P3E20.8)
 1320 FORMAT(' CF, EF, ZETAF =',1P3E20.8)
 1330 FORMAT(//28X,'EQUILIBRIUM CONSTANT'/,
     1 ' AKS, BKS, CKS, DKS, EKS =',1P5E18.8)
 1340 FORMAT(2(A9,2X,1PE12.5))
 1350 FORMAT(1X,A4,' STOICHIOMETRIC COEFS =',15I5/(26X,15I5))
 1360 FORMAT(1X,'3RD BODY EFFICIENCIES =',8F8.3/(29X,8F8.3))
 1370 FORMAT(1X,A4,'WARD SPECIES EXPONENTS =',8F8.3/(29X,8F8.3))
 1405 FORMAT(' FEE, FRE, TFLAGE =',1P3E20.8)
 1400 FORMAT(' CFE, EFE, ZETAFE =',1P3E20.8)
 1410 FORMAT(//11X,'FAST KINETIC AND EQUILIBRIUM REACTION',I5/,16X,A50/
     1 ' AS, BS, CS, DS, ES =',1P5E18.8)
 1420 FORMAT(' QEQ =',1PE18.8,10X,'NLM =',I5)
 1430 FORMAT(1X,A2,' =',15I5/(5X,15I5))
 1500 FORMAT(1X,'XINJ=',F10.5,' YINJ=',F10.5,' ZINJ=',F10.5)
 1510 FORMAT(2(6X,E18.8))
 1520 FORMAT(1X,'ANOZ=',E18.8,' T1INJ=',E18.8)
 1530 FORMAT(1X,'CONE=',F8.3,' TILTZ=',F8.3,' TILTX=',F8.3)
 1540 FORMAT(1X,'TSPMAS=',E18.8,' PMINJ=',E18.8)
 1550 FORMAT(1X,'AVERAGE INJECTION VELOCITY = ',E12.6,' CM/SEC',
     1     /,1X,'VDISP=',F8.3,' TPINJ=',F8.3)
 1560 FORMAT(6X,F10.3)
 1570 FORMAT(1X,F10.8,' CM TO ',F10.8,' CM --- ',E18.8)
 1600 FORMAT(1X,'ROPI=',E12.6,' TMMI=',E12.6,' TMLI=',E12.6,
     1      /1X,'LATENT HEAT=',E12.6)
 1605 FORMAT(1X,'PARTICLE SPECIES MASS FRACTIONS',5(2X,E12.6))
 1610 FORMAT(1X,'PSPHS=',E12.6,' PSPHL=',E12.6,' EMISSIVITY=',E12.6)
 1615 FORMAT(1X,'RPMIN=',E12.6,' RPMAX=',E12.6)
 1620 FORMAT(1X,'CARRIER GAS RATE OF SPECIES',I3,' =',5(2X,E12.6))
C==============================================================RINPUT
      END
